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Abstract: This paper presents and experiments approaches to solve a new 
bi-objective routing problem called the ring star problem. It consists of locating 
a simple cycle through a subset of nodes of a graph while optimizing two kinds 
of cost. The first objective is the minimization of a ring cost that is related to 
the length of the cycle. The second one is the minimization of an assignment 
cost from non-visited nodes to visited ones. In spite of its obvious bi-objective 
formulation, this problem has always been investigated in a single-objective way. 
To tackle the bi-objective ring star problem, we first investigate different stand- 
alone search methods. Then, we propose two cooperative strategies that com- 
bines two multiple objective metaheuristics: an elitist evolutionary algorithm 
and a population-based local search. We apply this new hybrid approaches 
to well-known benchmark test instances and demonstrate their effectiveness in 
comparison to non-hybrid algorithms and to state-of-the-art methods. 

Key-words: Multi-objective combinatorial optimization, evolutionary algo- 
rithms, local exploration, hybrid approaches, routing. 
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Metaheuristiques et leur hybridation pour 

resoudre le Ring Star Problem bi-objectif: 

une etude comparative 

Resume : Cc document presente et experimente differentes approches pour 
la resolution d'un nouveau probleme de routage bi-objectif appele le ring star 
problem. II consiste a localiser un cycle simple dans un sous-ensemble de 
noeuds d'un graphe en optimisant deux types de cout. Le premier objectif 
est la minimisation d'un cout d'anneau qui est lie a la longueur du cycle. Le 
deuxieme est la minimisation d'un coiit d'assigncmcnt depuis les noeuds non 
visites vers les noeuds visites. En depit de sa formulation bi-objectif evidente, 
ce probleme a toujours ete etudie de fagon mono-objectif. Pour attaquer le 
ring star problem bi-objectif, nous etudions d'abord differentes mcthodcs de re- 
cherche. Ensuite, nous proposons deux strategies cooperatives qui combinent 
deux metaheuristiques multi-objectif : un algorithme evolutionnaire elitiste et 
une recherche locale a base de population. Nous appliquons ccs approches hy- 
brides des instances de test reputees et demontrons leur efficacite par rapport 
aux algorithmes non-hybrides et a des methodes classiques de la litterature. 

Mots-cles : Optimisation combinatoirc multi-objectif, algorithmes cvolution- 
naires, exploration locale, approches hybrides, routage. 
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1 Introduction 

First introduced by Labbe et al. [3S], the Ring Star Problem (RSP), also known 
as Median Cycle Problem (MCP), is a recent routing problem. It aims to locate 
a simple cycle (the ring) containing a subset of nodes of a graph (including 
the depot), while nodes that do not belong to the ring must be connected to a 
visited one. Two costs are considered here, what naturally leads to a bi-criteria 
formulation. The first one is the ring cost associated to the ring itself and the 
second one is the assignment cost associated to the arcs directed from non- visited 
nodes to visited ones. Even so, the RSP has always been investigated in a single- 
objective way; either where both criterion are combined, or where a criteria is 
treated as a constraint. Thus, despite its numerous industrial interests, this is 
the first time that such a problem is regarded as a multi-objective one, probably 
due to its complexity. Indeed, it is highly combinatorial as, once is decided which 
nodes have to be visited or not, a traveling salesman problem is still to be solved. 

Recently, we studied the effectiveness of existing metaheuristics for solving 
the bi-objective RSP [HP. In this paper, we provide a deeper analysis of this 
work, and propose new stand-alone and hybrid approaches. As an initial step, 
we investigate four metaheuristics to approximate the set of efficient solutions 
for the problem under consideration. First, IBMOLS, a population-based local 
search recently proposed in [3], is designed with a variable neighborhood. Next, 
we present steady-state variations of two well-known multi-objective search 
methods, namely IBEA ^SB] and NSGA-II [TJ]. Last, we propose a new simple 
elitist evolutionary algorithm, SEEA, able to solve any kind of multi-objective 
combinatorial optimization problems. We compare all these metaheuristics to 
each other on state-of-the-art benchmarks test instances. As a second step, we 
propose two new cooperative strategies between the local search method and 
SEEA, based on a periodic and on an auto-adaptive scheme, that try to benefit 
of the advantages of each method it is compound of. The performance of these 
hybrid metaheuristics is shown in comparison to non- hybrid approaches. 

The paper is organized as follows. Section [2] is devoted to the bi-objective 
RSP. Then, four metaheuristics for multi-objective optimization are introduced 
in Section [31 A general presentation of these methods and their application 
to the RSP is followed by a parameter analysis and a comparative case study. 
In Section m we propose new cooperation schemes to solve multi-objective op- 
timization problems. We experiment the resulting search methods on the bi- 
objective RSP and we compare the obtained computational results to the pre- 
vious ones. At last, conclusions and perspectives are drawn in the last section. 

2 The Bi-objective Ring Star Problem 

In this section, we first present some basic concepts, notation and definitions 
related to multi-objective optimization. Next, we provide a formulation of the 
Ring Star Problem (RSP) as a bi-objective problem. Finally, we survey the 
literature relating to the RSP and discuss its industrial concerns. 
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2.1 Multi-objective Optimization 

A Multi-objective Optimization Problem (MOP) aims to optimize a set of n > 2 
objective functions fi, f2, ■ ■ ■ , fn simultaneously. Each objective function can be 
either minimized or maximized. Let X denote the set of feasible solutions in the 
decision space, and Z the set of feasible points in the objective space. A decision 
vector X = {xi, X2, . . . , Xk) is represented by a vector of k decision variables. In 
the case of a Multi- objective Combinatorial Optimization Problem (MCOP), a 
decision vector x G X has a finite set of possible values. To each decision vector 
X £ X is assigned exactly one objective vector z G Z on the basis of a vector 
function f : X ^ Z with z = f{x) — {fi{x), f2{x), . . . , fn{x))- We will assume, 
throughout the paper, that objective values are normalized. To achieve this, the 
minimum and the maximum value of each objective function are used in order 
to adaptively replace each objective function by its corresponding normalized 
function, so that its values lie in the interval [0, 1]. Without loss of generality, 
we assume that Z C 3?" and that all n objective functions have to be minimized. 
Therefore, a MOP can be formulated as follows: 

(MOP) = ( '™™' /(^) = (/i(^)'/2(a;),-.-,/n(a;)) ^. 

^ ' \ subject to x e X ^ ' 



Definition 1 An objective vector z G Z weakly dominates another objective 
vector z' G Z if and only ifVi £ {1, 2, . . . , n}, Zi < z[. 

Definition 2 An objective vector z G Z dominateq^ another objective vector 
z' G Z if and only ifii £ {1, 2, . . . , n}, Zi < z[ and 3j £ {1, 2, . . . , n} such as 
z, < z]. 

Definition 3 An objective vector z £ Z is non-dominated if and only if there 
does not exist another objective vector z' £ Z such that z' dominates z. 

A solution a; £ X is said to be efficient (or Pareto optimal, non- dominated) if its 
mapping in the objective space results in a non-dominated point. The set of all 
efficient solutions is the efficient (or Pareto optimal) set, denoted by Xe- The 
set of all non-dominated vectors is the non- dominated front (or the trade-off 
surface), denoted by Z^. A common approach in solving MOPs is to find or 
to approximate the minimal set of efficient solutions, i.e. a solution x £ Xe for 
each non-dominated vector z £ Zjv such as f{x) — z (in case there exists multi- 
ple solutions mapping to the same non-dominated point). But, generating the 
entire set of Pareto optimal solutions is usually infeasible due to the complexity 
of the underlying problem or to the large number of optima. Therefore, the 
overall goal is often to identify a good approximation of it. Population-based 
metaheuristics are commonly used to this end as they naturally find multiple and 
well-spread non-dominated solutions in a single simulation run. A reasonable 
basic introduction to multi-objective optimization can be found in [13J. But, 
the interested reader could refer to [H [H] for more details about evolutionary 
multi-objective optimization and to |17j for more details about multi-objective 
combinatorial optimization. 



^We will also say that a decision vector x £ X dominates a decision vector x' G X ii f(x) 
dominates f{x'). 
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2.2 Problem Definition 

The Ring Star Problem (RSP) can be described as follows. Let G = {V,E,A) 
be a complete mixed graph where V = {ui, W2, . . . ,««} is a set of vertexes, 
E = {[vi,Vj]\vi,Vj G V,i < j} is a set of edges, and A = {(vi,Vj)\vi,Vj e V} 
is a set of arcs. Vertex vi is the depot. To each edge [i^ijUj] we assign a 
non- negative ring cost cij , and to each arc {vi , Vj ) we assign a non-negative 
assignment cost dij. The RSP consists of locating a simple cycle through a 
subset of nodes V' <ZV (with vi G V') while (i) minimizing the sum of the ring 
costs related to all edges that belong to the cycle, and (ii) minimizing the sum 
of the assignment costs of arcs directed from every non-visited node to a visited 
one so that the associated cost is minimum. An example of solution is given in 
Figure [H where solid lines represent edges that belong to the ring and dashed 
lines represent arcs of the assignments. 



O--, 



o. 




.o 



--0 



Figure 1 : An example of solution for the ring star problem. 
The first objective is called the ring cost and is defined as: 



E 



Ci-jXi 



(2) 



[vi,v.i]<£E 



where Xij is a binary variable equal to 1 if and only if the edge [vi, Vj] belongs 
to the cycle. The second objective, the assignment cost, can be computed as 
follows: 



E 

ViGV\V 



min di. 



(3) 



Let us remark that these two objectives are comparable only if we assume that 
the ring cost and the assignment cost are commensurate one to another, what is 
rarely the case in practice. Furthermore, the fact of privileging a cost compared 
to the other is closely related to the decision-maker preferences. However, the 
RSP is an NP-hard combinatorial problem since the particular case of visiting 
the whole set of nodes is equivalent to a traditional traveling salesman problem. 

2.3 Related Works 

The RSP belongs to the class of location-allocation problems aiming at locating 
structures in a graph (see [32' for a review). It has initially been formulated 
by Labbe et al. |33j in two different ways. In the first formulation (denoted 
MCPl in [3S] and more often called 'ring star problem'), the sum of both ob- 
jectives is to be optimized. In the second formulation (MCP2 [33], usually 
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called 'median cycle problem'), the ring cost is to be minimized, while the as- 
signment cost is bounded by a prefixed value. Even if it was not explicitly 
noticed by the original authors, these two kinds of formulation are commonly 
employed to convert a multi-objective problem into a single-objective one by us- 
ing two well-known scalar approaches; respectively an aggregation method and 
an epsilon-constraint method^. The first formulation of the problem has been 
more widely studied in [34j . The authors used a branch-and-bound method and 
successfully solved TSPLIB and randomly generated instances involving up to 
200 nodes in less than two hours. In [S^, the same authors investigated the 
second formulation of the problem using a very close method. Finally, one or 
both versions of the problem have been heuristically solved by a variable neigh- 
borhood tabu search [H], an evolutionary algorithm [37j, a multi-start greedy 
add heuristic [47j . and a variable neighborhood tabu search hybridized with a 
greedy randomized adaptive search procedure [15] . 

A very similar problem, namely the capacitated m-ring star problem, has 
been introduced in [2]. It extends the first formulation by bounding the num- 
ber of nodes that can be assigned to a ring using a prefixed parameter (the 
ring capacity). This problem has been exactly solved by a branch-and-bound 
method and experiments have been performed on random and real-world in- 
stances. Note that the same problem has also been tackled in ^1] by a hybrid 
metaheuristic. Another variant of the RSP is the Steiner Ring Star Problem, 
where only optional nodes can be visited [SHj. In addition, Beasley and Nasci- 
mento j5] introduced a Single Vehicle Routing-Allocation Problem (SVRAP) in 
which a non-visited node can be either assigned to the cycle or isolated. In 
this case, a corresponding isolated cost has to be minimized as well. Vogt et 
al. [S3] proposed a tabu search to solve a single-objective problem resulting of 
this formulation. 

As reflected in the survey of Jozefowiez et al. |27| . an increasing number 
of multi-objective routing problems appeared in the literature in recent years. 
However, in spite of its numerous industrial applications (see subsection 12. 4p . 
the RSP has never been explicitly investigated in a multi-objective fashion. 
Nevertheless, Current and Schilling [TI] defined two multiple criteria variants 
of a very similar problem: the median tour problem and the maximal covering 
tour problem. In both, one criterion is the minimization of the total length 
of the tour, while another one is the maximization of the access to the tour 
for non-visited nodes. To tackle these problems, the authors used a kind of 
lexicographic method. This approach consists of assigning a priority value to 
each objective function, and in solving the problem one objective after another. 
Also, Dorner et al. |16j recently formulated a tri-objective optimization problem 
of tour planning for mobile health care facilities in Senegal, closely related to the 
SVRAP [6] introduced earlier. A mobile facility has to visit a subset of nodes. 
Non-visited nodes are then assigned to their closest tour stop or are regarded 
as unable to reach a tour stop (within a predefined maximum distance). The 
objectives are (i) to minimize the ratio between medical working time and total 
working time, (ii) to minimize the average distance to the nearest tour stops 
and (Hi) to maximize a coverage criterion. To do so, a Pareto ant colony 
optimization algorithm as well as two genetic algorithms (namely VEGA [35] 
and MOGA (20J) were designed to solve real-world instances. 



The reader is referred to |43j for more details about scalar methods for solving MOPs. 
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2.4 Industrial Interests 

The RSP has a wide range of industrial interests, especially in telecommunica- 
tions. Of course, practical applications may contain additional constraints. For 
instance, as noticed by Labbe et al. [Ml [55] , a ring-based network is designed 
to interconnect a set of hubs in the case of digital data services [51]. Some 
concentrators are installed on a subset of locations and are interconnected on a 
ring network (the Internet) while the remaining locations are assigned to this 
concentrators (the Intranet). Very close problems arise in rapid transit systems 
planning [33] or while designing optical networks [5]. Also, other kinds of ap- 
plications appears in the postal collection or delivery routes design, where the 
distance between a customer and a collection point have to be reasonable. For 
instance, post-box location while taking both the collection cost and the user 
inconvenience into account has been studied in [31j . Besides, other applications 
closely related to the RSP are the location of circular shaped transportation 
infrastructure (such as metro lines or motorways), the location of recyclable 
garbage collection bins, and school bus routing. Finally, the routing of essential 
health care services, already investigated in [H[TB] among other authors, consists 
of a mobile clinic servicing an area without being able to visit every population 
nodes. Then, unvisited ones have to reach the nearest tour stop by their own 
to be medically treated. 

3 Metaheuristics to Solve the Bi-objective Ring 
Star Problem 

Four metaheuristics are proposed to tackle the bi-objective RSP: a variable 
neighborhood iterative population-based Local Search (LS) and three Evolu- 
tionary Algorithms (EAs) . These algorithms are respectively steady-state vari- 
ations of IBMOLS [31 IBEA [56] and NSGA-II [M], as weU as a new search 
method, SEEA, introduced in this paper for the first time. IBMOLS and IBEA 
are both recent indicator-based metaheuristics, whereas NSGA-II is one of the 
most often used Pareto-based resolution approach. In this section, RSP-specific 
components are described after we presented the main characteristics of the LS 
and of the EAs. 

3.1 A Multi-objective Local Search 

Since they are easily adaptable to the multi-objective context, many of the 
search algorithms proposed to tackle MOPs are EAs. However, LS algorithms 
are known to be effective metaheuristics for solving real-world applications [S] 
[22]. Several neighborhood search methods to solve MOPs have been proposed 
in the literature. Most of them are based on a set of aggregations or scalariza- 
tions of the objective functions, what is the case, for instance, in [2T ] [23 ] [26 ] [51] . 
Pareto-based LS algorithms, like the one proposed by Paquete et al. [35], are 
more rare. An Indicator-Based Multi-Objective Local Search (IBMOLS for 
short) has recently been presented in ^ . IBMOLS is a generic population-based 
multi-objective LS dealing with a fixed population size. This allows to obtain 
a set of efficient solutions in a single run without specifying any mechanism to 
control the number of solutions during the search process. Moreover, IBMOLS 
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presents an alternative to aggregation- and Pareto-based multi-objective meta- 
heuristics. Indeed, as proposed in [5S], it is assumed that the optimization goal 
is given in terms of a binary quality indicator I |58j which can be regarded as 
an extension of the Pareto dominance relation. A value I{A, B) quantifies the 
difference in quality between two approximated efficient sets A and B. So, if R 
denotes a reference set (which can be the Pareto-optimal set Xe or any other 
set), the overall optimization goal can be formulated as: 

arg min^gj^ /(A, R) , (4) 

where J7 denotes the space of all efficient set approximations. As noted in '56j, 
R does not have to be known, it is just required in the formalization of the 
optimization goal. Since R is fixed, / actually represents a unary function that 
assigns a real number reflecting the quality of each approximation set according 
to the optimization goal. If / is dominance preserving |56j, I{A, R) is minimum 
for A = R. One of the main advantages of indicator-based optimization is that 
no specific diversity preservation mechanism is generally required, according to 
the indicator being used. 

The IBMOLS algorithm maintains a population P of size N . Then, it gen- 
erates the neighborhood of a solution contained in P until a good solution is 
found (i.e. one that is better than at least one solution of P in terms of the 
indicator being used). By iterating this simple principle to every solution of P, 
we obtain a local search step. The whole local search procedure stops when the 
archive of potentially efficient solutions has not received any new item during 
a complete local search step. Moreover, as local search methods are usually 
performed in an iterative way, a population re-initialization scheme has to be 
designed after each local search. As noticed in 3J, several strategies can be 
used within an Iterative IBMOLS (I-IBMOLS). Solutions can be re- initialized 
randomly, and crossover or random noise can be applied to solutions belonging 
to the efficient set approximation. The interested reader could refer to [5] for 
more details about IBMOLS and I-IBMOLS. 

A beneficial feature of this LS is the low number of parameters that are 
required. In addition to the population size, the binary quality indicator to be 
used and the population re-initialization strategy (between each local search) 
are the only two other problem-independent parameters. Indeed, several quality 
indicators can be used within IBMOLS. Some examples can be found in [24, 211 
[58] . One of them is the binary additive e-indicator (/e+) [Ml [58], inspired by the 
concept of e-dominance introduced by Laumanns et al. [36]. This indicator is 
particularly well-adapted to indicator-based search and seems to be efficient on 
different kinds of problems (see, for instance, [3[5B])- It is capable of obtaining 
both a well-converged and a well-diversified Pareto front approximation. I^^ 
computes the minimum value by which a solution x G AT can be translated 
in the objective space to weakly dominate another solution x' G X. For a 
minimization problem, it is defined as follows: 

Ie+{x,x')^ max (/j(a;) -/,(x')) . (5) 

ie{l,...,n} 

Furthermore, to evaluate the quality of a solution a; G A according to a whole 
population P and a binary quality indicator /, and then to compute the fitness 
value of a;, different approaches exist. As proposed in |56j . we will here consider 
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an additive technique that amphfies the influence of solutions mapping to dom- 
inating points over solutions mapping to dominated ones, what can be outlined 
as follows: 

Iix,P\{x})^ Y, -e-^(-*'-)/% (6) 

x*£P\{x} 

where k > is a scaling factor. However, the initial experiments were not 
satisfactory because the algorithm was not able to find the extreme points of 
the trade-off surface. As pointed out in [H], this is known to be one of the 
drawbacks of the e-dominance relation, apparently due to the high convexity 
of the front. Indeed, the authors illustrate that one of the limitation of the 
e-dominance is that the extreme points of the Pareto front are usually lost. To 
tackle this problem, we add a condition preventing the deletion of solutions 
corresponding to the extreme non-dominated vectors during the replacement 
step of IBMOLS. Besides, the population re-initialization scheme used between 
each local search procedure is based on random noise, such as in the basic 
simulated annealing algorithm [5]. This noise consists of multiple mutations 
applied to N different randomly chosen solutions contained in the archive of 
potentially efficient solutions. If the size of the archive is less than N, the 
population is filled with random solutions. 

3.2 Multi-objective Evolutionary Algorithms 

The first two multi-objective EAs designed to tackle the bi-objective RSP are 
variations of two state-of-the-art search methods, namely IBEA 56 and NSGA- 
II |14j . Some minor modifications have been carried out to improve the algo- 
rithms for the particular case of the addressed problem, for which a preliminary 
investigation revealed that the set of non-dominated points was, in general, very 
large. Finally, a generic approach to solve MOPs, called SEEA, is proposed in 
this paper for the first time and is presented in details. 

3.2.1 IBEA 

Introduced by Zitzler and Kiinzli [56j, the Indicator-Based Evolutionary Algo- 
rithm (IBEA) is, like IBMOLS, an indicator-based metaheuristic. The fitness 
assignment scheme of this EA is based on a pairwise comparison of solutions con- 
tained in a population by using a binary quality indicator. As within IBMOLS, 
no diversity preservation technique is required, according to the indicator being 
used. The selection scheme for reproduction is a binary tournament between 
randomly chosen individuals. The replacement strategy is an environmental one 
that consists of deleting, one-by-one, the worst individuals, and in updating the 
fitness values of the remaining solutions each time there is a deletion; this is con- 
tinued until the required population size is reached. Moreover, an archive stores 
solutions mapping to potentially non-dominated points in order to prevent their 
loss during the stochastic search process. However, in our case, and in contrast 
to the IBEA defined in [H], this archive is updated at each generation since 
the beginning of the EA, so that the output size is not necessarily less than or 
equal to the population size. Just like for the IBMOLS algorithm, the indicator 
used within IBEA is the additive e-indicator; and the same mechanism has been 
used to prevent the loss of the extreme points on the trade-off surface. More 
information about this algorithm can be found in I56j . 
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3.2.2 NSGA-II 

At each generation of NSGA-II (the Non-dominated Sorting Genetic Algorithm 
II introduced by Deb et al. in [15 )i the solutions contained in the population 
are ranked into several classes. Individuals mapping to vectors from the first 
front all belong to the best efficient set; individuals mapping to vectors from the 
second front all belong to the second best efficient set; and so on. Two values 
are computed for every solution of the population. The first one corresponds to 
the rank the corresponding solution belongs to, and represents the quality of the 
solution in terms of convergence. The second one, the crowding distance, consists 
of estimating the density of solutions surrounding a particular point of the 
objective space, and represents the quality of the solution in terms of diversity. 
A solution is said to be better than another if it has a best rank value, or in 
case of equality, if it has the best crowding distance. The selection strategy is 
a deterministic tournament between two random solutions. At the replacement 
step, only the best individuals survive, with respect to the population size. 
Likewise, an external population is added to the steady-state NSGA-II in order 
to store every potentially efficient solution found during the search. The reader 
is referred to [M] for more details about NSGA-II. 

3.2.3 SEEA 

If evaluating a solution in the objective space is not too much time consum- 
ing (what is the case for our problem), computing fitness values and diversity 
information are generally the most computationally expensive steps of a multi- 
objective EA. Based on this observation, we propose here a simple search method 
for which none of these phases is required. The resulting EA, called Simple Eli- 
tist Evolutionary Algorithm (SEEA for short), is detailed in Algorithm [TJ An 
archive of potentially efficient solutions is updated at each generation, and the 
individuals contained in the main population are generated by applying variation 
operators to randomly chosen archive members. The replacement step is a gen- 
erational one, i.e. the parent population is replaced by the offspring one. Note 
that the initial population can, for instance, be filled with random solutions. 
Then, as proposed in [57] among other authors, the archive is not only used as 
an external storage, but it is integrated into the optimization process during the 
selection phase of the EA, what is called elitism. The use of elitism is an im- 
portant issue within evolutionary multi-objective optimization [37' and SEEA 
is in somehow related to other elitist multi-objective EAs such as PAES 3D], 
PESA [lOj or SEAMO [52J. But, contrary to other approaches, no strategy 
to preserve diversity or to manage the size of the archive is involved here, as 
solutions are selected randomly and the archive is unbounded. The biggest ad- 
vantage of this EA is that the population (or the population size if solutions 
are randomly initialized) is the only problem-independent parameter. If non- 
dominated solutions are relatively close to each other in the decision space and 
if the archive is not too small compared to the main population, we believe that 
SEEA may convergence to a good approximation of the efficient set in a very 
short run-time. However, in some case, this method may prematurely converge 
or may appear inefficient if promising solutions are far from each other. 
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Algorithm 1 Simple Elitist Evolutionary Algorithm (SEEA) 



Input: P Initial population 

Output: A EfScicnt set approximation 



Step 1: 
Step 2: 

Step 3: 

Step 4 
Step 5. 
Step 6. 



Initialization. A <~ non-dominated individuals oi P; N ^- \P\; P' <— 0. 

Selection. Repeat until \P'\ = N: randomly select an individual from A and add 

it to the offspring population P' . 

Variation. Apply crossover and mutation operators to individuals of the offspring 

population P' . 

Replacement. P ^ P'; P' ^%. 

Elitism. A <— non-dominated individuals oi AVJ P. 

Termination. If a stopping criteria is satisfied return A, else go to Step 2. 



3.3 Application to the Ring Star Problem 

This section provides the problem-specific steps of the metaheuristics introduced 
earlier. Components designed for the particular case of the bi-objective RSP, 
such as the encoding mechanism, the population initialization as well as the 
neighborhood, mutation and crossover operators, are described below. 

3.3.1 Solution Encoding 

The encoding mechanism used to represent a RSP solution, for both the LS 
and the EAs, is based on the random keys concept proposed by Bean [5j. This 
implementation has already been successfully applied for a single-objective ver- 
sion of the RSP in [47] . To each node Vi belonging to the ring we assign exactly 
one random key xi £ [0, 1[, where xi = 0. A special value is assigned to un- 
visited nodes. Thus, the ring route associated to a solution corresponds to 
the nodes ordered according to their random keys in an increasing way; i.e. if 
Xi < Xj, then Vj comes after Vi. As an example, a possible representation for 
the cycle (vi,V7,t'4,U9,U2,W6) is given in Figure [21 vertexes W3,W5,U8 and viq are 
assigned to a visited node in such a way that the associated assignment cost is 
minimum. 



Vertex 


VI 


V2 


V3 


V4, 


V5 


V(i 


v-j 


vs, 


vg 


vio 


Random key 





0.7 


- 


0.3 


- 


0.8 


0.2 


- 


0.5 


- 



Figure 2: A RSP solution represented by random keys. 



3.3.2 Population Initialization 

For each optimization method, the initial population has been generated ran- 
domly. Each node has a probability p — 0.5 that it will be visited or not, and 
to each visited vertex we associate a key value randomly generated between 
and 1. 
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3.3.3 Neighborhood and Mutation Operators 

As the RSP is both a routing problem and an assignment problem, different 
neighborhood and mutation operators have to be designed. Here, we consider 
the following: 

• insert operator: adds an unvisited node vt in the cycle, the position where 
to insert Vi is chosen in order to minimize the ring cost 

• remove operator: removes a vertex Vj of the ring 

• 2-opt operator: applies a 2-opt operator between two nodes of the cycle Vi 
and Vj, i.e. the sequence of visited nodes located between Vi and Vj is 
reversed. 

For the LS, the neighbors of a solution are randomly explored, without consid- 
ering any order between these three operators; and each neighbor is at most 
visited once. Moreover, note that it is not necessary to completely re-evaluate 
a solution each time a neighborhood operator is applied. Thus, after an in- 
sert neighborhood operator, we just have to re-assign unvisited nodes in order 
to minimize the assignment cost. After a remove neighborhood operator, we 
just have to re-assign the nodes that were previously assigned to the one that 
has been removed. And, after a 2-opt neighborhood operator, we just have to 
recompute the ring cost, the assignment cost being unchanged. In the case of 
mutations, the operators are applied to randomly chosen vertexes. 

3.3.4 Crossover Operator 

The crossover operator is a quadratic crossover closely related to the one pro- 
posed in [37]. Two randomly selected solutions si and S2 are divided in a partic- 
ular position. Then, the first part of si is combined with the second part of S2 to 
build a first offspring, and the first part of S2 is combined with the second part 
of si to build a second offspring. Every node retains its random key so that 
it enables an easy reconstruction of the new individuals. Figure [3] illustrates 
a recombination between two solutions {vi,V'r,V4,vg,V2,VQ) and {vi,vs,V4,V3,v<^) 
after vertex vq, what gives rise to a couple of new solutions {vi,vs,V4,V2,vg) 
and (wi,W7,wg,W4,W3,W5). Thanks to the random keys encoding mechanism, so- 
lutions having a different ring size can easily be recombined, even if the initial 
ring structures are generally broken in the offspring solutions. 



Parent 1 



Parent 2 



Vertex 



Random key 
Vertex 



Random key 



vi 





vi 



V2 



0.7 

V2 



V3 



V4 



1)3 



0.8 



0.3 

V4 



V5 



V6 



0.7 



0.9 



0.8 

V6 



vr 



0.2 

V7 



Vg 



Vs 



0.2 



0.5 

vg 



t'lO 



I'lO 



Offspring 1 
Offspring 2 



Vertex 


VI 


V2 


V3 


V4 


«5 


V6 


vr 


vs 


vg 


fio 


Random key 





0.7 


- 


0.3 


- 


0.8 


- 


0.2 


- 


- 



Vertex 



Vl V2 vz Vi V. 



Random key - 0.8 0.7 0.9 

Figure 3: Crossover operator. 
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Another crossover operator preserving a bigger part of the initial ring struc- 
tures, has been experimented, but it was tending to reduce the number of nodes 
belonging to the cycle and then was causing a premature convergence with so- 
lutions having a small number of visited nodes. 



3.4 Experiments 

The metaheuristics described in the previous section have all been implemented 
using the ParadisEO-MOEO librarjll [39). ParadisEO-MOEO is a C++ white- 
box object-oriented framework dedicated to the reusable design of metaheuris- 
tics for multi-objective optimization. All the algorithms share the same base 
components for a fair comparison between them. Computational runs were 
performed on an Intel Core 2 Duo 6600 (2 x 2.40 GHz) machine, with 2 GB 
RAM. 



3.4.1 Experimental Protocol 

Benchmark Test Instances. Search methods have all been tested on differ- 
ent benchmark instances taken from the TSPLIEq |35] . These instances involve 
between 50 and 300 nodes. The number at the end of an instance's name repre- 
sents the number of nodes for the instance under consideration. Let Uj denote 
the distance between two nodes Vi and Vj of a TSPLIB file. Then, the ring 
cost Cij and the assignment cost dij have both been set to Uj for every pair of 
nodes Vi and Vj. Moreover, note that, for each resolution method proposed to 
tackle the bi-objective RSP, the search process stops after a certain amount of 
run time. As shown in Table [l] this stopping criteria has been arbitrary set 
according to the size of the instance under consideration. 



Table 1: Stopping criteria: running time. 



Instance 


Running 




time 


eil51 


20" 


st70 


1' 


kroAlOO 


2' 


bier 127 


5' 



Instance 


Running 




time 


kroAlSO 


10' 


kroA200 


20' 


pr264 


30' 


pr299 


50' 



Performance Assessment. For each one of the search method, a set of 20 
runs per instance, with different initial populations, has been performed. In 
order to evaluate the quality of the non-dominated front approximations ob- 
tained for a specific test instance, we follow the protocol given by Knowles et 
al. [H]. For a specific instance, we first compute a reference set Zj^ of non- 
dominated points extracted from the union of all these fronts. Second, we 
define 2;™°^ = (z^""^, z™""^), where z^°'^ (respectively z™"^) denotes the upper 
bound of the first (respectively second) objective in the whole non-dominated 
front approximations. Then, to measure the quality of an output set A in com- 
parison to Z*^, we compute the difference between these two sets by using the 



^ParadisEO is available at http://paradiseo.gforge.lnria.fr 

'jhttp : //www . iwr . uni-heidelberg ■ de/groups/comopt/sof tware/TSPLIB95/ 1 
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r 



■ Reference set Z,, 
• Approximation set S 



n 



"^ 



Figure 4: Illustration of the hypervolume difference {Ij^) between a reference 
set Z^ and a non-dominated front approximation S (shaded area). 



unary hypervolume metric [57], (1.05 x z™"^, 1.05 x z™"^) being the reference 
point. As illustrated in Figure SI the hypervolume difference indicator (I^) 
computes the portion of the objective space that is weakly dominated by Z^ 
and not by A. The more this measure is close to 0, the better is the approxima- 
tion A. Furthermore, we also consider the additive e- indicator proposed in [58j . 
Contrary to the one proposed in Equation [5l this indicator is used to compare 
non-dominated set approximations, and not solutions. The unary additive e- 
indicator {ll_^_) gives the minimum factor by which an approximation A has to 
be translated in the criterion space to weakly dominate the reference set Z'^. 
ll^ can be defined as follows: 



where 



/,V(A)=/,+ (A,Z]V), (7) 



I^+{A,B) = min{Vz e B,3z' e A : z'^ - e < z„\fl < i < n} . (8) 



As a consequence, for each test instance, we obtain 20 hypervolume differences 
and 20 Ie_|_ measures, corresponding to the 20 runs, per algorithm. As sug- 
gested by Knowles et al. [55^, once all these values are computed, we perform 
a statistical analysis on pairs of optimization methods for a comparison on a 
specific test instance. To this end, we use the Mann- Whitney statistical test 
as described in [W, with a p- value lower than 5%. Hence, for a specific test 
instance, and according to the p-value and to the metric under consideration, 
this statistical test reveals if the sample of approximation sets obtained by a 
search method Mi is significantly better than the sample of approximation sets 
obtained by a search method Af2, or if there is no significant difference between 
both. Note that all the performance assessment procedures have been achieved 
using the performance assessment tool suite provided in PlSi^lfl [?]■ 



^The package is available at |http : //www ■ t ik ■ ee . ethz . ch/pisa/assessment . html [ 
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3.4.2 Local Search Parameter Analysis 

Parameter Setting. For the IBMOLS algorithm, we set the population size 
according to the number of vertexes involved in the instance (see Table [5]) . For 
each instance, a small (S), a medium (M) and a large (L) population size have 
been experimented. Likewise, the noise rate for the population re-initialization 
in the iterated version of IBMOLS (denoted by I-IBMOLS) is set to a fixed 
percentage of the instance's size. We investigate three different values for this 
noise rate: 5%, 10% and 20%. Then, 0.05 x n, 0.1 x n and 0.2 x n random 
mutations are applied respectively for a problem with n nodes. Finally, the 
scaling factor k of the Equation [H] is set to 0.05. 



Table 2: Population size parameter for the I-IBMOLS algorithm. 



Instance 


Population 


size 




S 


M 


L 


eil51 


10 


20 


30 


stlO 


10 


20 


30 


kroAlOO 


20 


30 


50 


bierl27 


20 


30 


50 


kroAlSO 


20 


30 


50 


kroA200 


20 


30 


50 


pr264 


30 


50 


70 


pr299 


30 


50 


70 



Computational Results and Discussion. The results we obtained for the 
I^ and the ll^ metrics are respectively summarized in Table [3] and Table H) 
These tables give, for every instance, the number of times that a specific para- 
metrization is significantly worse than another one for the I-IBMOLS algorithm. 
So, a lower score is better, and zero means that the parametrization under 
consideration is never significantly outperformed. According to the I^ metric, 
a medium population size with a random noise of 5% or 10% seems to perform 
well on every test instance, except for the pr264 instance where a medium 
population size and a random noise of 10% is significantly worse than another 
parametrization. Moreover, as shown in Tabled a medium population size with 
a random noise parametrized to 10% generally seems to perform slightly better 
than with a random noise of 5% according to the I^^ metric. The only instances 
for which such a parametrization is significantly worse than at least another one 
are the kroA200 and the pr264 problems. Then, for all the other instances, we 
choose to set the number of individuals in the population to the medium size and 
the noise rate to 10% for comparing I-IBMOLS to the other metaheuristics. For 
the kroA200 instance, the same parametrization will be used as it seems to be 
the best compromise between the two metrics. Finally, for the pr264 instance, 
we will use a small population size and a noise rate of 20%. The value of all the 
parameters used by I-IBMOLS for the remainder of the paper are summarized 
in Table m 
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Table 3: Performance comparison of different population size and noise rate 
parameter settings for 1-IBMOLS according to the 1^ metric by using a Mann- 
Whitney statistical test with a significance level of 5%. Each column gives the 
number of algorithms with another parameter setting by which it is statistically 
outperformed for the instance under consideration. 



noise rate 




5% 






10% 






20% 




pop. size 


S 


M 


L 


s 


M 


L 


S 


M 


L 


eil51 


3 





3 


3 





3 


8 


1 


3 


stlO 


6 








7 








5 





1 


kroAlOO 


1 





6 


1 





6 


2 


2 


6 


bterl27 


1 





6 


1 





6 


5 


1 


6 


km A 150 


1 





6 


4 





6 


5 


1 


6 


kroA200 


2 





1 


2 





2 


8 


2 


2 


pr264 


5 





2 


2 


1 


4 





1 


5 


pr299 


6 








8 





3 


3 





4 



Table 4: Performance comparison of different population size and noise rate 
parameter settings for I-IBMOLS according to the I^_|_ metric by using a Mann- 
Whitney statistical test with a significance level of 5%. Each column gives the 
number of algorithms with another parameter setting by which it is statistically 
outperformed for the instance under consideration. 



noise rate 




5% 






10% 






20% 




pop. size 


S 


M 


L 


S 


M 


L 


S 


M 


L 


etlSl 


7 





3 


6 





3 


5 





3 


stlO 


7 








7 








6 





4 


kroAlOO 








6 








6 


1 


2 


6 


bierl27 








6 


1 





6 


1 


1 


6 


km A 150 


2 


1 


6 


3 





6 


2 





6 


kmA200 


8 


1 


1 


7 


1 


1 


1 





1 


pr264 


2 


2 


6 


1 


2 


6 





2 


6 


pr299 


8 


3 


3 


6 





3 








6 
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3.4.3 Evolutionary Algorithms Parameter Analysis 

Parameter Setting. For all the EAs, population sizes of 50, 100 and 200 
individuals have been experimented. Note that the instance-dependent popula- 
tion size setting used for I-IBMOLS has also been investigated, but the results 
were overall less efficient. The scaling factor k used by IBEA is set to 0.05. At 
last, the crossover probability is set to 0.25, and the mutation probability to 
1.00, with a probability of 0.25, 0.25 and 0.50 for the remove, the insert and 
the 2-opt operator, respectively. 

Computational Results and Discussion. As pointed out in Table [51 the 
results obtained by IBEA arc quite heterogeneous from a test instance to an- 
other, and no general conclusion can be drawn. Thus, we will use a population 
size of 50 for the eil51, the st70, the pr264 and the pr299 instances, a population 
size of 100 for the kroAlOO and the bier 127 instances, and a population size of 
200 for the kroAlSO and the kroA200 instances. For NSGA-II, the best perfor- 
mances are usually obtained with a large population size. Then, a population 
of 100 individuals for instances with less than 100 nodes and of 200 individuals 
for instances with more than 100 nodes will be used for the following. These pa- 
rameters are never significantly outperformed according to the metrics we used. 
Finally, as shown in Table [3 the population size does not appear to have a big 
influence on the outcome of SEEA, what is not surprising given the design of the 
algorithm (see Section [5. 2. 3p . As a consequence, we will use a population of 100 
individuals as this is the only parametrization that is strictly never significantly 
outperformed for every test instance. All the parameters used by the EAs for 
the next section are recapitulated in Table [S] 



Table 5: Performance comparison of different population size parameter settings 
for IBEA according to the I^ and the I^_,_ metrics by using a Mann- Whitney 
statistical test with a significance level of 5%. Each column gives the num- 
ber of algorithms with another parameter setting by which it is statistically 
outperformed for the instance under consideration. 







^H 






IcV 




pop. size 


50 


100 


200 


50 


100 


200 


eil51 








2 








2 


stlO 








2 





1 


1 


kroAlOO 








1 











bier 127 




















kroAlSO 


1 








1 








kroA200 





2 





1 


2 





pr264 





1 


1 





1 


1 


pr299 





1 


1 





1 


1 



3.4.4 Comparison between Search Methods 

Previous experiments allowed us to determine satisfying parameters for each 
one of the proposed algorithms. The parameter setting is summarized in Ta- 
ble [SI In this section, experiments were carried out with the aim of comparing 
the results obtained by the different metaheuristics to solve the problem under 



RR n° 6515 



Liefooghe et al. 



Table 6: Performance comparison of different population size parameter set- 
tings for NSGA-II according to the I^ and the ij^ metrics by using a Mann- 
Whitney statistical test with a significance level of 5%. Each column gives the 
number of algorithms with another parameter setting by which it is statistically 
outperformed for the instance under consideration. 







Ih 






I.\ 




pop. size 


50 


100 


200 


50 


100 


200 


eil51 


1 





2 


1 





2 


stlO 


2 








2 








kroAlOO 


2 







2 







bier 127 


2 







2 







km A 150 


2 







2 







kroA200 


2 







2 







pr264 


2 







2 







pr299 


2 







2 








Table 7: Performance comparison of different population size parameter settings 
for SEEA according to the I^ and the I^^ metrics by using a Mann- Whitney 
statistical test with a significance level of 5%. Each column gives the num- 
ber of algorithms with another parameter setting by which it is statistically 
outperformed for the instance under consideration. 
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Table 8: Parameter setting used to compare search methods. 





I-IBMOLS 


IBEA 


NSGA-II 


SEEA 


Instance 


pop. size 


noise rate 


pop. size 


pop. size 


pop. size 


eil51 


M = 20 


10% 


100 


100 


100 


st70 


M = 20 


10% 


50 


100 


100 


kroAlOO 


M = 30 


10% 


100 


200 


100 


bier 127 


M = 30 


10% 


100 


200 


100 


kroAlSO 


M = 30 


10% 


200 


200 


100 


kroA200 


M = 30 


10% 


200 


200 


100 


pr264 


S = 30 


20% 


50 


200 


100 


pr299 


M = 50 


10% 


50 


200 


100 
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consideration. A further analysis of these results was conducted in comparison 
to the previous ones. Then, Tableland Table fTUl present the results we respec- 
tively obtained for the I^ and the ij^ metrics. These tables show, for every test 
instance, the outcome of the statistical test as well as the significance level given 
by the p-value. With respect to the I^ metric, we can see that TIBMOLS is 
never outperformed by any other search method for all the instances we tested. 
The only instance for which I-IBMOLS does not perform better than the other 
algorithms is pr264, where there is no significant difference between the results 
it obtained and the ones of SEEA. Furthermore, according to the same metric, 
SEEA seems overall more efficient than IBEA and NSGA-II. Indeed, it always 
performs significantly better than NSGA-II and is statistically outperformed by 
IBEA only on the bierl27 and the kroA200 instances. Lastly, IBEA performs 
significantly better than NSGA-II on every test instance. Likewise, as shown in 
Table [101 the results we obtained for the I^^ metric, reveal that I-IBMOLS is 
globally more efficient than SEEA, IBEA and NSGA-II on most test instances. 
Nevertheless, it is not the case on the st70 and kroAlOO instances, where there 
is no significant difference between I-IBMOLS and SEEA, and on the pr264 in- 
stance, where I-IBMOLS is statistically outperformed by SEEA. Moreover, the 
results of NSGA-II are always worst than the ones of the other methods for all 
the instances. At last, SEEA provides better performance than IBEA on every 
problem, except for the kroA200 instance, where the difference between both 
is not significant, and for the bierl27 instance, where IBEA performs better. 
Then, to summarize the results given in Table El and Table Uni for the I^ and 
the Ig_|_ metric respectively, I-IBMOLS seems to be the overall most efficient 
optimization method to solve the bi-objective ring star problem. Then comes 
SEEA that performs better than any other EA. Finally, IBEA performs better 
than NSGA-II, the later being always outperformed by the other search meth- 
ods. All these observations are confirmed by the box-plots given in Figure[5]and 
FigureEl We can see that the measures obtained by NSGA-II and, to a lesser ex- 
tent, IBEA are both generally quite far from the ones obtained by IBMOLS and 
SEEA, except for some particular instances such as bier 127 and kroA200. More- 
over, the differences between the measures obtained by I-IBMOLS and SEEA 
are very thin for most of the cases, even if SEEA is generally outperformed, as 
pointed out above. 

In order to study the differences between I-IBMOLS and SEEA more finely, 
examples of the empirical attainment function obtained by these two algorithms 
are represented in Figure[5]and FigurefTUlfor the bier 127 and the pr264 instances 
respectively. Theses figures illustrate the limit of the objective space that is 
attained by at least 90% of the runs by both methods. The reader is referred 
to [19j for more details about the concept of empirical attainment function. 
For the bierl27 instance, where the results obtained by I-IBMOLS are much 
better than the ones of SEEA with respect to both metrics, we can see that 
SEEA seems to have more trouble in finding solutions having a low ring cost, 
whereas it seems slightly better than I-IBMOLS for finding interesting solutions 
with a small assignment cost. For the pr264 instance, where SEEA is overall 
more efficient, it is less obvious to observe any superiority graphically. Such a 
result mainly seems to be due to the discontinuity of the trade-off surface, for 
which a good approximation is given in Figure [5] An example of a more often 
encountered Pareto front is illustrated in Figure [7] for the &«eri ^7 instance. 
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Table 9: Performance comparison for I-IBMOLS, IBEA, NSGA-II and SEEA 
according to the I^ metric by using a Mann- Whitney statistical test. The "p- 
value" columns give the p- value of the statistical test. The "T" columns give the 
outcome of the statistical test for a significance level of 5%: either the results of 
the search method located at a specific row are significantly better than those of 
the search method located at a specific column (;^), either they are worse (-(), 
or there is no significant difference between both (=). 
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Table 10: Performance comparison for I-IBMOLS, IBEA, NSGA-II and SEEA 
according to the lj_,_ metric by using a Mann- Whitney statistical test. The "p- 
value" columns give the p- value of the statistical test. The "T" columns give the 
outcome of the statistical test for a significance level of 5%: either the results of 
the search method located at a specific row are significantly better than those of 
the search method located at a specific column (;^), either they are worse (-(), 
or there is no significant difference between both (=). 
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Figure 5: Performance comparison for I-IBMOLS, IBEA, NSGA-II and SEEA 
according to the I^ and the ij^ metric (1). 



INRIA 



Metaheuristics and Hybridization for the Bi-objective Ring Star Problem 23 



Instance 



kroAlSO 



o- 
d 










^ 


-r- 












• 




































<_j - 
d 


1 






I.IBMOLS IBEA 


NSGA.II 


SE'EA 



kroA200 









o - 








i- 

d 

o 
o - 
d 






ni 
















o 


— '— 



pr264 





















o 

O - 

o 
d 

o 

O - 

o 




o 


1 1 










o 

O - 

o 


1 1 









pr299 



o- 












o 








en 




o - 




o 


-^ 






;3; 






d 


1 1 







Figure 6: Performance comparison for I-IBMOLS, IBEA, NSGA-II and SEEA 
according to the I^ and the I^_|_ metric (2). 
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Figure 7: Best set of non-dominated points found for the bierl27 test instance. 
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Figure 8: 90%-attainnient surface plot obtained by the approximation sets found 
by I-IBMOLS and SEEA for the bier 127 test instance. 
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Figure 9: Best set of non-dominated points found for the pr264 test instance. 
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Figure 10: 90%-attainment surface plot obtained by the approximation sets 
found by I-IBMOLS and SEEA for the pr264 test instance. 
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One of the main characteristics of the problem under consideration seems 
to be the high number of points located in the trade-off surface, as shown in 
Figure [7] and Figure [H Then, after a certain number of iterations, a large 
part of the population involved in NSGA-II, IBEA and I-IBMOLS might map 
to potentially non-dominated points. This could explain the low efficiency of 
NSGA-II. Since the same rank is assigned to the major part of the population, 
only the crowding distance is used to compare solutions. However, it is not 
the case in SEEA because all non-dominated solutions contained in the archive 
can potentially take part in the evolution engine. Moreover, the indicator-based 
fitness assignment scheme is obviously much more suited to discriminate poten- 
tially efficient solutions than the single crowding distance. Moreover, the high 
performance of I-IBMOLS in comparison to IBEA might depends on how close 
are the solutions which map to non-dominated points in the decision space. If 
these solutions are close to each other, according to the neighborhood opera- 
tors, a LS is known to be particularly well-suited to find additional interesting 
solutions by exploring the neighborhood of a potentially efficient solution. On 
the contrary, an EA usually explores the decision space in a more random way. 
Thus, a landscape analysis could be interesting to study the bi-objective RSP 
in more depth. Besides, the results obtained by SEEA are quite satisfying in 
contrast to those obtained by the other EAs. Indeed, NSGA-II and IBEA are 
two state-of-art search methods for multi-objective optimization and are widely 
used in the community. Nevertheless, despite its simplicity, SEEA seems to be 
much more efficient for solving the B-RSP in a given amount of time. It could 
then be interesting to experiment the ability of this method for the resolution 
of other kinds of multi-objective combinatorial optimization problems. 

As a next step, we will try to solve problem instances involving an even 
bigger number of nodes to verify if our observations are still valid. It could then 
be interesting to design a cooperation scheme between two different methods 
(i.e. the local search procedure and an evolutionary algorithm) in order to 
benefit from the respective quality of each one of them. The resulting hybrid 
metaheuristic could be particularly effective for solving large size problems. 

4 A Cooperative Approach 

This section presents a cooperative approach combining SEEA and the non- 
iterative version of IBMOLS in order to solve the bi-objective RSP. Two variants 
are proposed: a periodic one that operates a systematic cooperation and an 
auto-adaptive one that decides on-line when the cooperation must occur. 

4.1 Motivation 

Designing metaheuristics for solving combinatorial optimization problems is 
generally a matter of intensification and diversification. This is even more pro- 
nounced for MOPs where the goal is to find a well-converged and well-diversified 
efficient set approximation. However, LSs are known to be particularly efficient 
as intensifying methods whereas EAs are clearly powerful to explore the deci- 
sion space thanks to their variation operators. Instead of trying to improve one 
method in term of diversification or the other in term of intensification, a com- 
mon approach is to hybridize both in order to make them cooperate and then to 
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benefit of their respective behaviors. Thus, hybrid metaheuristics have shown 
their efficiency to solve different kinds of optimization problems |49j , including 
MOPs [H] . They generally consist of a search method that cooperates with a 
second one in order to improve a solution or a set of solutions. 

In the previous section, we saw that SEEA and I-IBMOLS were overall more 
efficient than NSGA-II and IBEA to approximate the efficient set for the prob- 
lem under consideration. Moreover, these two methods are quite different to 
each other and do not explore the search space in the same way. SEEA has 
been conceived in order to find a rough approximation of the Pareto set in a 
very short amount of time whereas the non-iterative version of IBMOLS is able 
to improve an approximated set in a very efficient way. It could then be in- 
teresting to design a cooperation scheme between these two algorithms. The 
resulting hybrid metaheuristic could be particularly efficient for solving large 
size problems. Furthermore, both methods maintain a secondary population 
(the archive) in parallel of the main population to store non-dominated solu- 
tions. This archive is not only used as an external storage (what is the case, for 
instance, in IBEA and NSGA-II), but also takes part in the evolution engine as 
it serves to build new solutions to explore. Thus, each method can manage its 
own population and therefore use the archive as a single shared memory. 

4.2 Cooperative Schemes 

The general idea of our hybridization scheme is to run SEEA and to lunch IB- 
MOLS regularly by using a subset of archive items as an initial population. As 
the non-iterative version of IBMOLS stops when its own archive does not re- 
ceive any new efficient solution anymore, we can restart the SEEA process until 
the next step of the hybrid algorithm. A step of the hybrid metaheuristic can, 
for instance, be defined by a certain amount of time or by a certain number of 
generations. Besides, as SEEA uses the non-dominated solutions found by IB- 
MOLS to create new ones and vice versa, the global archive is the only memory 
shared by the two search agents to exchange information. Resulting from this, 
we can imagine two versions of the hybrid algorithm: (i) a periodic version, 
in which IBMOLS is launched at each step, and {ii) an auto-adaptive version, 
in which IBMOLS is launched at a specific step only if a condition is verified. 
These two approaches will be denoted by PCS (for Periodic Cooperative Search) 
and ACS (for Auto-adaptive Cooperative Search) in the remainder of the paper, 
and are respectively illustrated in Figure [Tl] and Figure [TH 

The ACS method decides by itself, and on-line, if it would be interesting 
to launch IBMOLS at a specific instant of the search process. The condition 
that we consider here to start IBMOLS during a particular step of the ACS is 
that the archive of potentially efficient solutions does not improve enough re- 
garding to the optimization scenario. One possibility is to measure the quality 
of the current archive A^ in comparison to the one of the previous step A^~^. 
Different metrics exist to evaluate the convergence of an approximated efficient 
set regarding to another. For instance, we could have computed the hyper- 
volume [ST] of A* and A^^^ and measure the difference between both. But, 
even if it evaluates the quality of an approximated efficient set both in term of 
convergence and diversity, the hypervolume metric has the drawback of being 
time consuming. Thus, let us introduce the contribution metric C proposed 
by Meunier et al. [42] . This metric gives an idea of the quality of an approxi- 
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Cooperative Search (PCS). adaptive Cooperative Search (ACS). 



mated efficient set S\ compared to another one 5*2 in term of convergence and 
can be computed in a reasonable running time. If 5** denotes the set of non- 
dominated solutions given by 5*1 U 5*2, the contribution of 5*1 on 5*2 evaluates 
the proportion of non-dominated solutions represented by 5*1 in S* . Then, if 
/(5'i) = /(S'2), C(5i,S'2) = C(S'2,S'i) = 0.5. If every solution of S^ is domi- 
nated by at least one solution of S'2, CiSi^S-i) — 0. And, generally speaking, 
C{S\,S-2) -I- C(S'2,5'i) = 1. In our case, at each step t of the ACS, we com- 
pute the contribution of the current archive A* on the archive of the previous 
step A*~^. Thus, as non-dominated solutions are not lost between two steps, A* 
is at least as good as A*+^ and the set of non-dominated solutions of A* U A*~^ 
is then A*. As a consequence, we know that C(A*,^*~^) e [0.5, 1]. Assuming 
that the archive does not improve enough if the contribution of yl* on yl*^^ is 
less than a pre-defined threshold 5 6 [0.5, 1], we choose to launch IBMOLS only 
if C(j4*, A*~^) < 5. Let us remark that an ACS with ^8 — 0.5 would result in 
SEEA, and that an ACS with a 5 = 1.0 would result in PCS, the time spent in 
calculating the different contribution values in less. 

4.3 Related Works 

Different schemes exist on how two search methods can be combined. According 
to the taxonomy proposed in 49J , two levels (low and high) and two modes (relay 
and cooperative) of hybridization can be distinguished. In low-level hybridiza- 
tion, a given function of a metaheuristic is replaced by another metaheuristic. 
In high-level hybridization, there is no direct relationship between the inter- 
nal workings of a metaheuristic. A second hierarchical classification deals with 
the mode of cooperation. The relay mode consists of a method applied after 
another one in a pipeline way, the last one using the output of the previous 
one as its input. On the contrary, teamwork hybridization represents a coop- 
erative optimization model in which search agents exchange information with 
the others. Here, we consider that the optimization process alternates between 
two search algorithms: SEEA and the non-iterative version of IBMOLS. Both 
methods are self-contained and cooperates with each other via a global archive. 
Therefore, the resulting hybrid metaheuristic can be classified on the High-level 
Teamwork Hybrid (HTH) class of the taxonomy introduced in '49] and can then 
be denoted by IITII(SEEA-I-IBMOLS). In their survey on hybrid metaheuris- 
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tics to solve multi-objective combinatorial optimization problems [18j . Ehrgott 
and Gandibleux identify three categories of methods hybridizing an EA with a 
neighborhood search algorithm: (i) an hybridization to make a method more 
aggressive, (ii) an hybridization to drive a method and (Hi) an hybridization 
for exploiting complementary strengths. The last one consists of alternating 
between both search methods, what is the case within our hybridization. But, 
most of existing approaches occur in a pipeline way, like in the relay mode. 
The few teamwork hybridization techniques of the literature to solve MOPs are 
most often compound of the same kind of metaheuristics, what is generally the 
case in the island model. Moreover, this class of hybrid methods often uses a 
neighborhood search method combining a set of scalarizing functions. The orig- 
inality of the approach proposed in this paper is that search agents are based 
on different types of metaheuristics, are hybridized in a teamwork mode, and do 
not use any scalar approach to convert the multiple objective functions into a 
single one. Furthermore, one of the variant we propose, the auto-adaptive one, 
automatically detects when to start the neighborhood search according to the 
optimization scenario. 

4.4 Experiments 

In order to experiment the efhciency of our two cooperative approaches, we 
compare them to SEEA and to the iterative version of IBMOLS by using the 
same experimental protocol as the one defined in Section 13.4.11 After having 
set the parameters for each hybrid search method, we give some computational 
results and we discuss the contribution of the hybridization for the resolution of 
the B-RSP. Finally, we give an idea of the behavior of our approach compared 
to a state-of-art single-objective approach proposed in [M] . 

4.4.1 Parameter Setting 

For both cooperative search methods, the population size managed by SEEA is 
set to 100, and the population size managed by IBMOLS is set according to the 
instance under consideration. These sizes have been set on the same way that we 
did for the stand-alone iterative version of IBMOLS on the previous section (see 
Table [2]). But, for large-size instances, initial experiments were not satisfying 
since the hybrid algorithms were generally not able to launch IBMOLS more 
than once during the search process as it was too much time consuming. For 
this reason, we bounded the IBMOLS population size to 30. Then, a IBMOLS 
population of 20 individuals has been set for instances with less than 100 nodes, 
and a IBMOLS population of 30 individuals has been set for instances with 100 
nodes and more. The step t of the hybrid algorithms has been set to 0.5% of 
the total run time for the instance under consideration, rounded to the nearest 
upper integer. Finally, we investigated different S values for the ACS method: 
0.6, 0.7, 0.8 and 0.9. As pointed out above, an ACS with a S value of 0.5 and 
1.0 is similar to PCS and SEEA, respectively. 

4.4.2 Computational Results and Discussion 

ACS Parameter Analysis. As a first step, we analyze the influence of the 
6 parameter on the results of the ACS method. Results are summarized in Ta- 
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ble[TT] A value of 0.6 is generally less efficient than at least one other setting on 
most test instances according to both metrics. For the other values, the results 
are roughly comparable, even if a (5 = 0.8 is never statistically outperformed by 
any other value according to I^ and I^, . Therefore, we will use this parameter 
setting for the next experiments. 

Table 11: Performance comparison of different S parameter settings for the 
ACS method according to the I^ and the I^_|_ metrics by using a Mann- Whitney 
statistical test with a significance level of 5%. Each column gives the num- 
ber of algorithms with another parameter setting by which it is statistically 
outperformed for the instance under consideration. 
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Comparison between Search Methods. In addition to the benchmark test 
instances we investigated previously, we experimented two bigger ones: pr439 
and prl002. The maximum running time for these new instances has been set 
to 70 and 100 minutes, respectively. Following the general trend identified in 
Section I3.4.2[ we set the instance-specific parameters for the iterative version of 
IBMOLS as follows: a population of 70 and 100 individuals respectively, and 
a noise rate of 10%. The other parameters were set in the same way than for 
other instances. 

Table [T^ and Table [T2] respectively summarize the results we obtained by 
comparing I-IBMOLS, SEEA, PCS and ACS according to the 1;^^ and to the 
Ij^ metrics. For the 10 benchmark test instances we experimented, we can see 
that both cooperative search methods are never statistically outperformed by 
SEEA. Moreover, they perform significantly better according to at least one 
metric on all instances except eil51. Similarly, PCS and ACS are at least as 
good as I-IBMOLS on every test instance. While there is generally no significant 
difference for small-size instances, the efficient set approximations found by the 
cooperative approaches are statistically better than I-IBMOLS for problems 
with 150 nodes and more (except for the pr299 instance). Thus, contrary to 
what has been pointed out in the previous section about SEEA, the cooperative 
methods do not seem to have trouble in finding solution having a low ring 
cost in comparison to I-IBMOLS, as illustrated in the empirical attainment 
function given in Figure [Tsl At last, in most cases, the differences between PCS 
and ACS arc not statistically significant, except for the bierl27 instance where 
PCS performs better, and for the kroAlSO and the pr439 instances where ACS 
performs better according to the I^ metric. Box-plots for I^ and I^^ are given 
in Figure [T3] and Figure 1141 In general, these figures graphically confirm the 
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superiority of the cooperative search methods for test instances having a high 
number of nodes. Thus, to summarize the results obtained by the hybrid search 
methods, the benefit of the cooperation scheme is being felt in most large-size 
test instances, but the addition of an auto-adaptive aspect does not seem to have 
a major influence on the results. Nevertheless, as pointed out in Table [HI the 
difference between the average number of times that IBMOLS is launched during 
the search process of PCS and ACS is relatively thin. This can be explained by 
the fact that («) IBMOLS puts more time to find non-dominated solutions by 
starting with a population of poorer quality, what is the case in PCS comparing 
to ACS (at least for the first launch), and {ii) a part of the run-time allocated 
to the algorithm is used to compute a contribution value at every step of ACS, 
whereas PCS devotes all of its run-time in the search process. These two aspects 
lead to the fact that the number of times that IBMOLS is launched is, in the 
end, more or less balanced between both cooperative methods. 

Comparison with Exact Single-objective Results. As a last step, we 
provide a comparison between the results found for the bi-objective version 
of the RSP proposed in this paper and the ones of the mono-objective RSP 
investigated in [;34 , where both costs are summed up. Note that it was not 
possible to compare our results to the ones of the other formulation of a single- 
objective RSP investigated in [3511131 137], where the assignment cost is subject 
to a constraint, due to the way the bound has been fixed. In 33], the authors 
propose an exact method to solve a single-objective RSP aiming at minimizing 
the sum of the ring cost and of the assignment cost. But, in order to provide 
optimal solutions visiting approximately 25, 50, 75 and 100% of the total number 
of nodes, they set the ring cost Cij and the assignment cost dij between two 
nodes Vi and Vj in the following way: c^ — \o:lij~\ and dij = [(10 — o:)lij~\ with 
a G {3,5,7,9}, where lij denotes the distance between Vi and Vj given in the 
TSPLIB files. In order to make our results comparable to their, we transform 
the objective vectors of the best found Pareto front Z^ we obtained for every 
instance into a scalar value as detailed above, and we compare the minimal one 
with the optimunij found in [34 . Additionally, we compare the best solution 
visiting every nodes we have found with the optimal TSP solution available in 
the TSPLIB websitql|. Table \TE\ gives the error ratio between the best known 
value and the best one we have found for every identified (single) objective and 
every benchmark test instance. In comparison to |34j . this ratio is always under 
1.5%, and is mostly below 1% for every instance. The optimum is even found 
for the kroAlOO and the kroAlSO instances with an a = 9, and a better solution 
is found for the kroA200 instances with an a = 3 and an a = 9. As regards to 
the optimal TSP solutions (corresponding to solutions visiting every node in the 
case of the RSP), our results are similar for instances with 200 nodes or less. For 
larger instances, they are slightly less good, especially for the prl002 instance, 
where the error ratio is close to 10%. But, proportionally to the number of 
nodes, the running time available to solve this test instance is lower than for 
the other instances. To sum up, compared with single-objective optimal or near 
optimal results, the search methods we proposed in this paper to solve the B- 



®In fact, the authors imposed a time limit for their experiments. They report the best 
solution found so far for the instances exceeding this time limit, what is the case for the 
kroA200 instance with an o = 3, 5 and 9 in Table [TSl 

'|http : //www . iwr . uni-heidelberg ■ de/groups/comopt/sof tware/TSPLIBQsT] 
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Table 12: Performance comparison for I-IBMOLS, SEEA, PCS and ACS ac- 
cording to the IJi metric by using a Mann- Whitney statistical test. The "p- 
value" columns give the p- value of the statistical test. The "T" columns give 
the outcome of the statistical test for a significance level of 5%: either the re- 
sults of the search method located at a specific row are significantly better than 
those of the search method located at a specific column (;^), either they are 
worse (-<), or there is no significant difference between both (=). 
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Table 13: Performance comparison for I-IBMOLS, SEEA, PCS and ACS ac- 
cording to the I^_|_ metric by using a Mann- Whitney statistical test. The "p- 
value" columns give the p- value of the statistical test. The "T" columns give the 
outcome of the statistical test for a significance level of 5%: either the results of 
the search method located at a specific row are significantly better than those of 
the search method located at a specific column (>-), cither they are worse (-(), 
or there is no significant difference between both (=). 
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Figure 13: Performance comparison for I-IBMOLS, SEEA, PCS and ACS ac- 
cording to the I^ and the I^_,_ metric (1). 
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Figure 14: Performance comparison for I-IBMOLS, SEEA, PCS and ACS ac- 
cording to the I^ and the I^_^ metric (2). 
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Table 14: Average number of times that IBMOLS is launched during the search 
process of PCS and ACS. 
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Figure 15: 90%-attainnicnt surface plot obtained by the approximation sets 
found by PIBMOLS and PCS for the bierl27 test instance. 
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RSP are quite promising with regard to the relatively small computational time 
available and to the size of the problem instances to be solved. Nevertheless, 
it is true that the solutions used to compare the results are the best ones we 
obtained during the set of all experiments we performed, and may not be as 
good for a single simulation run of a given search method. 

Table 15: Error ratio between the ring cost value of the best known TSP solution 
and the best single objective value found in |34] in comparison to the best value 
found during our experiments. 
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5 Conclusion and Perspectives 

A new multi-objective routing problem, the bi-objective ring star problem, has 
been investigated in this paper for the first time. It aims to locate a cycle 
through a subset of nodes of a graph while minimizing a ring cost, related to 
the length of the cycle, and an assignment cost, from non- visited nodes to visited 
ones. Despite it is clearly bi-objective, this problem has always been tackled in a 
single-objective formulation, either where both costs are combined (34] . or where 
one cost is regarded as a constraint [351117]. As a first step, four population- 
based metaheuristics have been proposed to approximate the minimal set of 
efficient solutions for the problem under consideration. A population-based local 
search method, IBMOLS, recently proposed in [3], has first been designed, in an 
iterative way, with a variable neighborhood. Then, two state-of-the-art multi- 
objective evolutionary algorithms, namely IBEA [56J and NSGA-II pjlj, as well 
as a third one, the simple elitist evolutionary algorithm (SEEA), proposed here 
for the first time, have all been fitted for the particular case of the bi-objective 
ring star problem to be solved. Experiments were conducted using a set of 
multiple benchmark test instances. After having studied the influence of some 
parameters on the efficiency of the different search methods, we compared them 
to each other and we concluded that the iterative version of IBMOLS was, in 
general, significantly better than all the evolutionary algorithms on the problem 
instances we tackled. Nevertheless, the behavior of a simple search method like 
SEEA in comparison to the ones of IBEA and NSGA-II is very encouraging 
with regard to combinatorial problem solving. As a second step, we designed 
two cooperative schemes, able to solve any kind of multi-objective optimization 
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problems, between SEEA and the non-iterative version of IBMOLS. Both consist 
of using SEEA as the main search process and of launching IBMOLS regularly 
to intensify the search at a given time. The first hybrid metaheuristic, the 
periodic cooperative search (PCS), operates a systematic launching of IBMOLS 
at each of its step. The second one, the auto-adaptive cooperative search (ACS), 
launches IBMOLS at a given step of the algorithm only if a certain condition 
is verified. Thus, ACS evolves adaptively according to the search scenario and 
decides by itself, and online, when the cooperation must occur. In comparison 
to stand-alone metaheuristics, these two hybrid search methods statistically 
improves the results on a large number of test instances, and particularly on 
large-size ones. However, the efficiency difference between PCS and ACS is 
almost negligible. This can be explained by the fact that ACS spends time to 
compute whether the cooperation should occur at a given step of the algorithm. 
This is not the case for PCS, so that the latter devotes all of its computational 
time to the search process, even if a cooperation may start at an inopportune 
time, and then be less gainful. 

Although the approaches proposed in this paper are already promising, even 
compared to the mono- objective approaches of the literature, a few future re- 
search directions are open. First is the possibility to improve the population 
initialization strategy used within all the proposed search methods. Indeed, the 
simple one that is used in the paper is a bit rudimentary, as each initial solution 
has approximately half of its nodes that belong to the cycle. Surely, an initial 
population having a well-spread number of visited nodes among its members 
should accelerate the convergence of the search methods toward a well-spread 
set of non-dominated solutions. Moreover, the use of a traveling salesman prob- 
lem heuristic would be helpful in order to improve the ring cost of these initial 
solutions. Second, we pointed out that the crossover operator designed for the 
problem under consideration has a tendency to break the parent ring structures 
in the offspring individuals. Then, as proposed in [17], we could employ the 
same heuristic to improve the ring cost of the newly generated solutions. Third, 
given that several mutation operators are used within the evolutionary algo- 
rithms, it is always difficult to determine the appropriate probability selection 
of each one of them. Then, a possibility would be to set these rates adaptively, 
during the search process, according to the efficiency of each operator. Such a 
strategy has already been proposed in [4] in the frame of multi-objective opti- 
mization. Note that the last two points could greatly improve the efficiency of 
the evolutionary algorithms compared to the local search method, what would 
also be beneficial for the hybrid methods. From a purely algorithmic point of 
view, it would be interesting to try other quality indicators than the e-indicator 
within the fitness assignment scheme of IBMOLS. We know that this binary in- 
dicator, although experiencing some difficulties concerning the extreme points 
of the trade-off surface, performs in general much better than others |3j. But 
the indicator-based fitness assignment strategy using the unary hypervolume in- 
dicator recently proposed in j55| could be an interesting alternative to explore. 
Next, we saw that the stand-alone metaheuristics taking part in the cooperative 
approach proposed in this paper, SEEA and IBMOLS, can both be launched 
independently, exchanging information via the archive containing the set of non- 
dominated solutions found so far. Then, analogously to what has been proposed 
in COSEARCH [50i for the mono-objective case, these two search agents can be 
launched in parallel while an adaptive memory stores different kind of informa- 
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tion, including the archive. But, a similar generic model adapted to the case of 
multi-objective optimization still has to be established. Finally, we found out 
that SEEA was a good alternative to state-of-the-art evolutionary algorithms 
for multi-objective combinatorial optimization when a given, relatively small, 
amount of computational time is allowed. In the future, we plan to tackle other 
kinds of combinatorial problems within SEEA to verify if our observations are 
still valid, especially in the case where more than two objectives are involved. 
The same remark can be done with regard to the cooperative approaches we 
proposed. 
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